/****************************************************************************
 * Copyright (C) 2009-2010 SciTouch LLC
 * 
 * This file is part of Indigo toolkit.
 * 
 * This file may be distributed and/or modified under the terms of the
 * GNU General Public License version 3 as published by the Free Software
 * Foundation and appearing in the file LICENSE.GPL included in the
 * packaging of this file.
 * 
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 ***************************************************************************/

#include "base_cpp/scanner.h"
#include "base_cpp/tlscont.h"
#include "base_cpp/auto_ptr.h"

#include "molecule/molfile_loader.h"
#include "molecule/molecule.h"
#include "molecule/molecule_stereocenters.h"
#include "molecule/molecule_3d_constraints.h"
#include "graph/graph_highlighting.h"

MolfileLoader::MolfileLoader (Scanner &scanner) : 
_scanner(scanner),
TL_CP_GET(_stereo_care_atoms),
TL_CP_GET(_stereocenter_types),
TL_CP_GET(_stereocenter_groups),
TL_CP_GET(_bond_directions),
TL_CP_GET(_ignore_cistrans)
{
   reaction_atom_mapping = 0;
   reaction_atom_inversion = 0;
   reaction_atom_exact_change = 0;
   reaction_bond_reacting_center = 0;
   highlighting = 0;
   _rgfile = false;
   ignore_stereocenter_errors = false;
   treat_x_as_pseudoatom = false;
}

void MolfileLoader::loadMolecule (Molecule &mol, bool query)
{
   mol.clear();
   _readHeader(mol);

   if (_rgfile && !query)
      throw Error("RGfiles are allowed only for queries");

   _readCtabHeader();

   if (query)
   {
      mol.createQueryAtoms();
      mol.createQueryBonds();
   }

   if (_v2000)
   {
      _readCtab2000(mol, query);

      if (_rgfile)
         _readRGroups2000(mol);
   }
   else
   {
      _readCtab3000(mol, query);
      
      if (query)
         _readRGroups3000(mol);
   }

   _postLoad(mol, query);
}

void MolfileLoader::loadCtab3000 (Molecule &mol, bool query)
{
   _readCtab3000(mol, query);
   _postLoad(mol, query);
}

void MolfileLoader::_readHeader (Molecule &mol)
{
   if (_scanner.lookNext() == '$')
   {
      _rgfile = true;        // It's RGfile
      _scanner.skipString(); // Skip $MDL REV  1   Date/Time
      _scanner.skipString(); // Skip $MOL
      _scanner.skipString(); // Skip $HDR
   }

   // Skip header
   _scanner.readString(mol.name, true);
   _scanner.skipString();
   _scanner.skipString();

   if (_rgfile)
   {
      _scanner.skipString(); // Skip $END HDR
      _scanner.skipString(); // Skip $CTAB
   }
}

void MolfileLoader::_readCtabHeader ()
{
   QS_DEF(Array<char>, str);

   _scanner.readString(str, false);

   BufferScanner strscan(str);

   _atoms_num = strscan.readIntFix(3);
   _bonds_num = strscan.readIntFix(3);

   try
   {
      char version[6];
      int chiral_int;
      
      strscan.skip(6);
      chiral_int = strscan.readIntFix(3);
      strscan.skip(19);
      strscan.read(5, version);
      strscan.skipString();

      version[5] = 0;

      if (strcmp(version, "V2000") == 0)
         _v2000 = true;
      else if (strcmp(version, "V3000") == 0)
         _v2000 = false;
      else
         throw Error("bad molfile version : %s", version);

      _chiral = (chiral_int != 0);
   }
   catch (Scanner::Error &)
   {
      _chiral = false;
      _v2000 = true;
   }
}

void MolfileLoader::_readCtab2000 (Molecule &mol, bool query)
{
   MoleculeRGroups *rgroups = 0;
   MoleculeRGroupFragment *rgroup_fragment = 0;

   _init();

   int k;

   QS_DEF(Array<char>, str);

   // read atoms
   for (k = 0; k < _atoms_num; k++)
   {
      // read coordinates
      float x = _scanner.readFloatFix(10);
      float y = _scanner.readFloatFix(10);
      float z = _scanner.readFloatFix(10);

      _scanner.skip(1);

      char buf[4] = {0, 0, 0, 0};
      int label = 0;
      int query_atom_type = 0;
      int isotope = 0;
      int charge = 0;
      bool have_pseudo_atom = false;

      // read atom label and mass difference
      _scanner.readCharsFix(3, buf);
      
      isotope = _scanner.readIntFix(2);

      if (buf[0] == 'R' && (buf[1] == '#' || buf[1] == ' '))
      {
         if (!query)
            throw Error("RGroups are allowed only for queries");

         if (mol.getRGroups() == 0)
            rgroups = mol.createRGroups();

         query_atom_type = QUERY_ATOM_RGROUP;
      }
      else if (buf[0] == 'A' && buf[1] == ' ')
      {
         label = -QUERY_ATOM_A;
         // will later become QUERY_ATOM_A or superatom
      }
      else if (buf[0] == 'X' && buf[1] == ' ' && !treat_x_as_pseudoatom)
      {
         if (!query)
            throw Error("'X' label is allowed only for queries");
         query_atom_type = QUERY_ATOM_X;
      }
      else if (buf[0] == 'Q' && buf[1] == ' ')
      {
         if (!query)
            throw Error("'Q' label is allowed only for queries");
         query_atom_type = QUERY_ATOM_Q;
      }
      else if (buf[0] == 'L' && buf[1] == ' ')
      {
         if (!query)
            throw Error("atom lists are allowed only for queries");
         query_atom_type = QUERY_ATOM_LIST;
      }
      else if (buf[0] == 'D' && buf[1] == ' ')
      {
         label = ELEM_H;
         isotope = 2;
      }
      else if (buf[0] == 'T' && buf[1] == ' ')
      {
         label = ELEM_H;
         isotope = 3;
      }
      else
      {
         for (int i = 0; i < 3; i++)
         {
            if (isspace(buf[i]))
            {
               buf[i] = 0;
               break;
            }
            if (isalpha(buf[i]))
               buf[i] = (i == 0) ? toupper(buf[i]) : tolower(buf[i]);
         }

         if (*buf == 0)
            throw Error("empty atom #%d", k);

         label = Element::fromString2(buf);

         if (label == -1)
            have_pseudo_atom = true;

         if (isotope != 0)
            isotope = Element::getDefaultIsotope(label) + isotope;
      }

      int idx = mol.addAtom(label);
      Atom &atom = mol.getAtom2(idx);

      if (have_pseudo_atom)
         mol.setPseudoAtom(idx, buf);

      atom.pos.set(x, y, z);
      atom.isotope = isotope;

      // read atom charge
      charge = _scanner.readIntFix(3);
      _assignCharge(charge, atom);

      _scanner.readString(str, false);

      BufferScanner rest(str);

      try
      {
         // skip atom stereo parity
         rest.skip(3);

         int stereo_care, hydrogens_count;

         hydrogens_count = rest.readIntFix(3);
         stereo_care = rest.readIntFix(3);

         if (query)
         {
            mol.getQueryAtom2(idx).type = query_atom_type;
            mol.getQueryAtom2(idx).hydrogens = hydrogens_count;
         }
         else
         {  
            if (hydrogens_count > 0)
               throw Error("only a query can have nonzero atom hydrogens count");
            if (query_atom_type > 0)
               throw Error("only a query can have query atom");
         }

         if (stereo_care)
            _stereo_care_atoms[idx] = 1;

         // read valence
         int valence = rest.readIntFix(3);

         if (valence > 0 && valence < 14)
         {
            atom.explicit_valence = true;
            atom.valence = valence;
         }

         // skip "HO designator" and 2 unused fields
         rest.skip(9);

         int aam = 0, irflag = 0, ecflag = 0;

         // read and set reaction parameters
         aam = rest.readIntFix(3);    // atom-to-atom mapping number
         irflag = rest.readIntFix(3); // inversion/retension flag,
         ecflag = rest.readIntFix(3); // exact change flag

         if (reaction_atom_mapping != 0)
            reaction_atom_mapping->at(idx) = aam;
         if (reaction_atom_inversion != 0)
            reaction_atom_inversion->at(idx) = irflag;
         if (reaction_atom_exact_change != 0)
            reaction_atom_exact_change->at(idx) = ecflag;
      }
      catch (Scanner::Error &)
      {
      }
   }

   int bond_idx;

   for (bond_idx = 0; bond_idx < _bonds_num; bond_idx++)
   {
      int beg, end, type, stereo, topology, query_type = 0;

      beg = _scanner.readIntFix(3);
      end = _scanner.readIntFix(3);
      type = _scanner.readIntFix(3);
      stereo = _scanner.readIntFix(3);

      int stereocenter_type = 0;

      if (stereo == 1 || stereo == 6)
      {
         if (_chiral)
            stereocenter_type = MoleculeStereocenters::ATOM_ABS;
         else
            stereocenter_type = MoleculeStereocenters::ATOM_AND;

         if (stereo == 1)
            _bond_directions[bond_idx] = MoleculeStereocenters::BOND_UP;
         else // stereo == 6
            _bond_directions[bond_idx] = MoleculeStereocenters::BOND_DOWN;
      }
      else if (stereo == 4)
      {
         _bond_directions[bond_idx] = MoleculeStereocenters::BOND_EITHER;
         stereocenter_type = MoleculeStereocenters::ATOM_ANY;
      }
      else if (stereo == 3)
         _ignore_cistrans[bond_idx] = 1;
      else if (stereo != 0)
         throw Error("unknown number for bond stereo: %d", stereo);

      if (stereocenter_type > 0 && _stereocenter_types[beg - 1] == 0)
      {
         _stereocenter_types[beg - 1] = stereocenter_type;
         _stereocenter_groups[beg - 1] = 1;
      }

      switch (type)
      {
         case BOND_SINGLE: 
         case BOND_DOUBLE:
         case BOND_TRIPLE:
         case BOND_AROMATIC:
            break;
         case QUERY_BOND_SINGLE_OR_DOUBLE:
            if (!query)
               throw Error("'single or double' bonds are allowed only for queries");
            query_type = type;
            type = 0;
            break;
         case QUERY_BOND_SINGLE_OR_AROMATIC:
            if (!query)
               throw Error("'single or aromatic' bonds are allowed only for queries");
            query_type = type;
            type = 0;
            break;
         case QUERY_BOND_DOUBLE_OR_AROMATIC:
            if (!query)
               throw Error("'double or aromatic' bonds are allowed only for queries");
            query_type = type;
            type = 0;
            break;
         case QUERY_BOND_ANY:
            if (!query)
               throw Error("'any' bonds are allowed only for queries");
            query_type = type;
            type = 0;
            break;
         default:
            throw Error("unknown bond type: %d", type);
      }

      int idx = mol.addBond(beg - 1, end - 1, type);
      QueryBond *query_bond = 0;

      if (query)
      {
         query_bond = &mol.getQueryBond2(idx);
         query_bond->type = query_type;
      }

      if (type == BOND_AROMATIC)
         mol.setBondAromatic(idx);

      _scanner.readString(str, false);

      BufferScanner rest(str);

      try
      {
         rest.skip(3); // not used

         // bond topology
         topology = rest.readIntFix(3);

         if (topology != 0)
         {
            if (query_bond != 0)
               query_bond->topology = topology;
            else
               throw Error("bond topology is allowed only for queries");
         }

         // read reaction parameters
         int reaction_info;

         reaction_info = rest.readIntFix(3);

         if (reaction_bond_reacting_center != 0)
            reaction_bond_reacting_center->at(idx) = reaction_info;
      }
      catch (Scanner::Error &)
      {
      }
   }

   if (highlighting != 0)
      highlighting->init(mol);

   // Special MDL rule of notation without chiral flag
   //if (n_stereocenters > 1 && chiral)
   //{
   //   for (atom_idx = 0; atom_idx < _atoms.size(); atom_idx++)
   //      if (_atoms[atom_idx].stereocenter == ATOM_STEREO_AND)
   //         _atoms[atom_idx].stereocenter = ATOM_STEREO_NCF;
   //}

   int n_3d_features = -1;

   // read groups
   while (!_scanner.isEOF())
   {
      char c = _scanner.readChar();

      if (c == 'G')
      {
         _scanner.skipString();
         _scanner.skipString();
         continue;
      }
      if (c == 'M')
      {
         _scanner.skip(2);
         char chars[4];

         chars[3] = 0;
         _scanner.readCharsFix(3, chars);

         if (strncmp(chars, "END", 3) == 0)
         {
            _scanner.skipString();
            break;
         }
         // atom list
         else if (strncmp(chars, "ALS", 3) == 0)
         {
            if (!query)
               throw Error("atom lists are allowed only for queries");

            int atom_idx;
            int list_size;
            char excl_char;

            _scanner.skip(1);
            atom_idx = _scanner.readIntFix(3);
            list_size = _scanner.readIntFix(3);
            _scanner.skip(1);
            excl_char = _scanner.readChar();
            _scanner.skip(1);

            atom_idx--;

            QueryAtom &query_atom = mol.getQueryAtom2(atom_idx);

            query_atom.type = (excl_char == 'F') ? QUERY_ATOM_LIST : QUERY_ATOM_NOTLIST;
            mol.getAtom2(atom_idx).label = 0;

            _scanner.readString(str, false);

            BufferScanner rest(str);

            int i;
            for (i = 0; i < list_size; i++)
            {
               int j;

               chars[0] = chars[1] = chars[2] = chars[3] = ' ';

               for (j = 0; j < 4; j++)
               {
                  // can not read 4 characters at once because
                  // sqlplus cuts the trailing spaces
                  if (!rest.isEOF())
                     chars[j] = rest.readChar();
                  else
                     break;
               }

               if (j < 1)
                  throw Error("atom list: can not read element #%d", i);

               if (chars[0] == 'A' && chars[1] == ' ')
               {
                  if (list_size == 1)
                     query_atom.type = QUERY_ATOM_A;
                  else
                     throw Error("'A' inside atom list, if present, must be single");
               }
               else if (chars[0] == 'Q' && chars[1] == ' ')
               {
                  if (list_size == 1)
                     query_atom.type = QUERY_ATOM_Q;
                  else
                     throw Error("'Q' inside atom list, if present, must be single");
               }
               else
               {
                  BufferScanner elemscan(chars, sizeof(chars));
                  query_atom.addToList(Element::read(elemscan));
               }
            }
         }
         // atom charge
         else if (strncmp(chars, "CHG", 3) == 0)
         {
            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               mol.getAtom2(atom_idx).charge = _scanner.readIntFix(3);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "RAD", 3) == 0)
         {
            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               mol.getAtom2(atom_idx).radical = _scanner.readIntFix(3);
            }
            _scanner.skipString();
         }
         // atom isotope
         else if (strncmp(chars, "ISO", 3) == 0)
         {
            int n = _scanner.readIntFix(3);

            while (--n >= 0)
            {
               _scanner.skip(1);

               int atom_idx = _scanner.readIntFix(3) - 1;

               _scanner.skip(1);

               mol.getAtom2(atom_idx).isotope = _scanner.readIntFix(3);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "SUB", 3) == 0)
         {
            if (!query)
               throw Error("substitution counts are allowed only for queries");

            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               mol.getQueryAtom2(atom_idx).substitution_count = _scanner.readIntFix(3);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "RBC", 3) == 0)
         {
            if (!query)
               throw Error("ring bond count is allowed only for queries");

            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               mol.getQueryAtom2(atom_idx).ring_bond_count = _scanner.readIntFix(3);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "UNS", 3) == 0)
         {
            if (!query)
               throw Error("unaturated atoms are allowed only for queries");

            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               int unsaturation = _scanner.readIntFix(3);

               mol.getQueryAtom2(atom_idx).unsaturation = (unsaturation != 0);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "$3D", 3) == 0)
         {
            if (!query)
               throw Error("3D features are allowed only for queries");

            if (n_3d_features == -1)
            {
               mol.create3dConstraints();
               n_3d_features = _scanner.readIntFix(3);
               _scanner.skipString();
            }
            else
            {
               n_3d_features--;

               if (n_3d_features < 0)
                  throw Error("3D feature unexpected");

               _read3dFeature2000(mol);
            }
         }
         else if (strncmp(chars, "AAL", 3) == 0)
         {
            if (!query)
               throw Error("rgroups are allowed only for queries");

            if (mol.getRGroups() == 0)
               rgroups = mol.createRGroups();

            int site_idx = _scanner.readIntFix(3) - 1;
            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               int att_type = _scanner.readIntFix(3);

               rgroups->setAttachmentOrder(site_idx, att_type - 1, atom_idx);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "RGP", 3) == 0)
         {
            if (!query)
               throw Error("rgroups are allowed only for queries");

            if (mol.getRGroups() == 0)
               rgroups = mol.createRGroups();

            int n = _scanner.readIntFix(3);

            while (n-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               int rgroup = _scanner.readIntFix(3) - 1;

               rgroups->addRGroupOnSite(atom_idx, rgroup);
            }
            _scanner.skipString();
         }
         else if (strncmp(chars, "LOG", 3) == 0)
         {
            if (!query)
               throw Error("rgroups logic are allowed only for queries");

            if (mol.getRGroups() == 0)
               rgroups = mol.createRGroups();

            // skip something
            _scanner.skip(3);

            _scanner.skip(1);
            int rgroup_idx = _scanner.readIntFix(3) - 1;
            _scanner.skip(1);
            int if_then = _scanner.readIntFix(3) - 1;
            _scanner.skip(1);
            int rest_h = _scanner.readIntFix(3);
            _scanner.skip(1);

            QS_DEF(Array<char>, occurence_str);

            RGroup &rgroup = rgroups->getRGroup(rgroup_idx);

            rgroup.if_then = if_then;
            rgroup.rest_h = rest_h;

            _scanner.readString(occurence_str, true);
            _readRGroupOccurenceRanges(occurence_str.ptr(), rgroup.occurence);
         }
         else if (strncmp(chars, "APO", 3) == 0)
         {
            if (!query)
               throw Error("rgroups attachment points are allowed only for queries");

            if (mol.getRGroupFragment() == 0)
               rgroup_fragment = mol.createRGroupFragment();

            int list_length = _scanner.readIntFix(3);

            while (list_length-- > 0)
            {
               _scanner.skip(1);
               int atom_idx = _scanner.readIntFix(3) - 1;
               _scanner.skip(1);
               int att_type = _scanner.readIntFix(3);

               if (att_type == -1)
                  att_type = 3;

               for (int att_idx = 0; (1 << att_idx) <= att_type; att_idx++)
                  if (att_type & (1 << att_idx))
                     rgroup_fragment->addAttachmentPoint(att_idx, atom_idx);
            }

            _scanner.skipString();
         }
         else
            _scanner.skipString();
      }
      else if (c == 'A')
      {
         char chars[4];
         QS_DEF(Array<char>, pseudo);

         _scanner.skip(2);
         chars[3] = 0;
         int atom_idx = _scanner.readIntFix(3);

         atom_idx--;
         _scanner.skipString();
         _scanner.readString(pseudo, true);
         _preparePseudoAtomLabel(pseudo);
         mol.setPseudoAtom(atom_idx, pseudo.ptr());
      }
      else
         _scanner.skipString();
      //  throw Error("unsupported molfile header: %s", chars);
   }

   for (k = 0; k < _atoms_num; k++)
   {
      if (mol.getAtom(k).label == -QUERY_ATOM_A)
      {
         if (!query)
            throw Error("'any' atoms are allowed only for queries");

         mol.getAtom2(k).label = 0;
         mol.getQueryAtom2(k).type = QUERY_ATOM_A;
      }
   }
}

void MolfileLoader::_assignCharge (int value, Atom &atom)
{
   switch (value)
   {
      case 1: atom.charge = 3; break;
      case 2: atom.charge = 2; break;
      case 3: atom.charge = 1; break;
      case 4: atom.radical = 2; break;
      case 5: atom.charge = -1; break;
      case 6: atom.charge = -2; break;
      case 7: atom.charge = -3; break;
   }
}

void MolfileLoader::_read3dFeature2000 (Molecule &mol)
{
   // read 3D feature ID (see MDL ctfile documentation)
   int feature_id = _scanner.readIntFix(3);

   _scanner.skipString();

   Molecule3dConstraints *constraints = mol.get3dConstraints();

   switch (feature_id)
   {
      case -1: // point defined by 2 points and distance
      {
         AutoPtr<Molecule3dConstraints::PointByDistance> constr;

         constr.create();

         _scanner.skip(6);
         constr->beg_id = _scanner.readIntFix(3) - 1;
         constr->end_id = _scanner.readIntFix(3) - 1;
         constr->distance = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -2: // point defined by 2 points and percentage
      {
         AutoPtr<Molecule3dConstraints::PointByPercentage> constr;

         constr.create();

         _scanner.skip(6);
         constr->beg_id = _scanner.readIntFix(3) - 1;
         constr->end_id = _scanner.readIntFix(3) - 1;
         constr->percentage = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -3: // point defined by point, normal line, and distance
      {
         AutoPtr<Molecule3dConstraints::PointByNormale> constr;

         constr.create();

         _scanner.skip(6);
         constr->org_id = _scanner.readIntFix(3) - 1;
         constr->norm_id = _scanner.readIntFix(3) - 1;
         constr->distance = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -4: // line defined by 2 or more points (best fit line if more than 2 points)
      {
         AutoPtr<Molecule3dConstraints::BestFitLine> constr;

         constr.create();

         _scanner.skip(6);

         int amount = _scanner.readIntFix(3);
         if (amount < 2)
            throw Error("invalid points amount in M $3D-4 feature");

         constr->max_deviation = _scanner.readFloatFix(10);
         _scanner.skipString();
         _scanner.skip(6);

         while (amount-- > 0)
            constr->point_ids.push(_scanner.readIntFix(3) - 1);

         _scanner.skipString();
         constraints->add(constr.release());
         break;
      }
      case -5: // plane defined by 3 or more points (best fit line if more than 3 points)
      {
         AutoPtr<Molecule3dConstraints::BestFitPlane> constr;

         constr.create();

         _scanner.skip(6);

         int amount = _scanner.readIntFix(3);

         if (amount < 3)
            throw Error("invalid points amount in M $3D-5 feature");

         constr->max_deviation = _scanner.readFloatFix(10);
         _scanner.skipString();
         _scanner.skip(6);

         while (amount-- > 0)
            constr->point_ids.push(_scanner.readIntFix(3) - 1);

         _scanner.skipString();
         constraints->add(constr.release());
         break;
      }
      case -6: // plane defined by point and line
      {
         AutoPtr<Molecule3dConstraints::PlaneByPoint> constr;

         constr.create();

         _scanner.skip(6);
         constr->point_id = _scanner.readIntFix(3) - 1;
         constr->line_id = _scanner.readIntFix(3) - 1;
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -7: // centroid defined by points
      {
         AutoPtr<Molecule3dConstraints::Centroid> constr;

         constr.create();
         _scanner.skip(6);

         int amount = _scanner.readIntFix(3);

         if (amount < 1)
            throw Error("invalid amount of points for centroid: %d", amount);

         _scanner.skipString();
         _scanner.skip(6);

         while (amount-- > 0)
            constr->point_ids.push(_scanner.readIntFix(3) - 1);

         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -8: // normal line defined by point and plane
      {
         AutoPtr<Molecule3dConstraints::Normale> constr;

         constr.create();
         _scanner.skip(6);
         constr->point_id = _scanner.readIntFix(3) - 1;
         constr->plane_id = _scanner.readIntFix(3) - 1;

         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -9: // distance defined by 2 points and range
      {
         AutoPtr<Molecule3dConstraints::DistanceByPoints> constr;

         constr.create();
         _scanner.skip(6);
         constr->beg_id = _scanner.readIntFix(3) - 1;
         constr->end_id = _scanner.readIntFix(3) - 1;
         constr->bottom = _scanner.readFloatFix(10);
         constr->top    = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -10: // distance defined by point, line and range
      {
         AutoPtr<Molecule3dConstraints::DistanceByLine> constr;

         constr.create();

         _scanner.skip(6);
         constr->point_id = _scanner.readIntFix(3) - 1;
         constr->line_id = _scanner.readIntFix(3) - 1;
         constr->bottom = _scanner.readFloatFix(10);
         constr->top = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -11: // distance defined by point, plane and range
      {
         AutoPtr<Molecule3dConstraints::DistanceByPlane> constr;

         constr.create();
         _scanner.skip(6);
         constr->point_id = _scanner.readIntFix(3) - 1;
         constr->plane_id = _scanner.readIntFix(3) - 1;
         constr->bottom = _scanner.readFloatFix(10);
         constr->top = _scanner.readFloatFix(10);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -12: // angle defined by 3 points and range
      {
         AutoPtr<Molecule3dConstraints::AngleByPoints> constr;

         constr.create();
         _scanner.skip(6);
         constr->point1_id = _scanner.readIntFix(3) - 1;
         constr->point2_id = _scanner.readIntFix(3) - 1;
         constr->point3_id = _scanner.readIntFix(3) - 1;
         constr->bottom = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         constr->top    = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -13: // angle defined by 2 lines and range
      {
         AutoPtr<Molecule3dConstraints::AngleByLines> constr;

         constr.create();
         _scanner.skip(6);
         constr->line1_id = _scanner.readIntFix(3) - 1;
         constr->line2_id = _scanner.readIntFix(3) - 1;
         constr->bottom = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         constr->top    = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -14: // angles defined by 2 planes and range
      {
         AutoPtr<Molecule3dConstraints::AngleByPlanes> constr;

         constr.create();
         _scanner.skip(6);
         constr->plane1_id = _scanner.readIntFix(3) - 1;
         constr->plane2_id = _scanner.readIntFix(3) - 1;
         constr->bottom = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         constr->top    = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -15: // dihedral angle defined by 4 points
      {
         AutoPtr<Molecule3dConstraints::AngleDihedral> constr;

         constr.create();
         _scanner.skip(6);
         constr->point1_id = _scanner.readIntFix(3) - 1;
         constr->point2_id = _scanner.readIntFix(3) - 1;
         constr->point3_id = _scanner.readIntFix(3) - 1;
         constr->point4_id = _scanner.readIntFix(3) - 1;
         constr->bottom = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         constr->top    = (float)(_scanner.readFloatFix(10) * M_PI / 180);
         _scanner.skipString();

         constraints->add(constr.release());
         break;
      }
      case -16: // exclusion sphere defines by points and distance
      {
         AutoPtr<Molecule3dConstraints::ExclusionSphere> constr;

         int allowed_atoms_amount;
         Array<int> allowed_atoms;

         constr.create();
         _scanner.skip(6);
         constr->center_id = _scanner.readIntFix(3) - 1;
         constr->allow_unconnected = (_scanner.readIntFix(3) != 0);
         allowed_atoms_amount = _scanner.readIntFix(3);
         constr->radius = (float)(_scanner.readFloatFix(10));

         if (allowed_atoms_amount > 0)
         {
            _scanner.skipString();
            _scanner.skip(6);

            while (allowed_atoms_amount-- > 0)
               constr->allowed_atoms.push(_scanner.readIntFix(3) - 1);
         }

         _scanner.skipString();
         constraints->add(constr.release());
         break;
      }
      case -17: // fixed atoms
      {
         _scanner.skip(6);
         int amount = _scanner.readIntFix(3);
         _scanner.skipString();
         _scanner.skip(6);

         while (amount-- > 0)
            mol.fixed_atoms.push(_scanner.readIntFix(3) - 1);

         _scanner.skipString();
         break;
      }
      default:
         throw Error("unknown 3D feature in createFromMolfile: %d", feature_id);
   }
}

void MolfileLoader::_readRGroupOccurenceRanges (const char *str, Array<int> &ranges)
{
   int beg = -1, end = -1;
   int add_beg = 0, add_end = 0;

   while (*str != 0)
   {
      if (*str == '>')
      {
         end = 0xFFFF;
         add_beg = 1;
      } else if (*str == '<')
      {
         beg = 0;
         add_end = -1;
      } else if (isdigit(*str))
      {
         sscanf(str, "%d", beg == -1 ? &beg : &end);
         while (isdigit(*str))
            str++;
         continue;
      } else if (*str == ',')
      {
         if (end == -1)
            end = beg;
         else
            beg += add_beg, end += add_end;
         ranges.push((beg << 16) | end);
         beg = end = -1;
         add_beg = add_end = 0;
      }
      str++;
   }
   if (end == -1)
      end = beg;
   else
      beg += add_beg, end += add_end;
   ranges.push((beg << 16) | end);
}

int MolfileLoader::_asc_cmp_cb (const int &v1, const int &v2, const void *context)
{
   return v2 - v1;
}

void MolfileLoader::_postLoad (Molecule &mol, bool query)
{
   MoleculeRGroups *rgroups = mol.getRGroups();

   int i;

   if (query)
   {
      for (i = mol.vertexBegin(); i < mol.vertexEnd(); i = mol.vertexNext(i))
      {
         QueryAtom &query_atom = mol.getQueryAtom2(i);

         // Update attachment orders for rgroup bonds if they were not specified explicitly
         if (query_atom.type == QUERY_ATOM_RGROUP && mol.getRGroups()->isRGroupAtom(i))
         {
            const Vertex &vertex = mol.getVertex(i);

            if (vertex.degree() == 1 && rgroups->getAttachmentOrder(i, 0) == -1)
            {
               rgroups->setAttachmentOrder(i, 0, vertex.neiVertex(vertex.neiBegin()));
            } else if (vertex.degree() == 2 &&
               (rgroups->getAttachmentOrder(i, 0) == -1 || rgroups->getAttachmentOrder(i, 1) == -1))
            {
               int nei_idx_1 = vertex.neiVertex(vertex.neiBegin());
               int nei_idx_2 = vertex.neiVertex(vertex.neiNext(vertex.neiBegin()));

               rgroups->setAttachmentOrder(i, 0, __min(nei_idx_1, nei_idx_2));
               rgroups->setAttachmentOrder(i, 1, __max(nei_idx_1, nei_idx_2));
            } else if (vertex.degree() > 2)
            {
               int j;

               for (j = 0; j < vertex.degree(); j++)
                  if (rgroups->getAttachmentOrder(i, j) == -1)
                  {
                     QS_DEF(Array<int>, nei_indices);

                     nei_indices.clear_resize(vertex.degree());

                     for (int nei_idx = vertex.neiBegin(); nei_idx < vertex.neiEnd(); nei_idx = vertex.neiNext(nei_idx))
                        nei_indices.push(vertex.neiVertex(nei_idx));

                     nei_indices.qsort(_asc_cmp_cb, 0);

                     for (int order = 0; order < vertex.degree(); order++)
                        rgroups->setAttachmentOrder(i, order, nei_indices[order]);
                     break;
                  }
            }
         }
      }
   }

   for (i = mol.edgeBegin(); i < mol.edgeEnd(); i = mol.edgeNext(i))
   {
      const Edge &edge = mol.getEdge(i);

      if (_stereo_care_atoms[edge.beg] && _stereo_care_atoms[edge.end])
      {
         if (!query)
            throw Error("stereo care box is allowed only for queries");

         if (_ignore_cistrans[i])
             continue;

         // in fragments such as C-C=C-C=C-C, the single bond
         // have both ends 'stereo care', but should not be considered
         // as 'stereo care' itself
         if (MoleculeCisTrans::isGeomStereoBond(mol, i, 0, true))
            mol.getQueryBond2(i).stereo_care = 1;
      }
   }

   _createStereocenters(mol);

   mol.cis_trans.build(mol, _ignore_cistrans.ptr());

   if (query && mol.getRGroups() != 0)
   {
      int n_rgroups = rgroups->getRGroupCount();
      for (i = 0; i < n_rgroups; i++)
         if (rgroups->getRGroup(i).occurence.size() == 0 && rgroups->getRGroup(i).fragments.size() > 0)
            rgroups->getRGroup(i).occurence.push((1 << 16) | 0xFFFF);
   }

   mol.have_xyz = true;
}

void MolfileLoader::_readRGroups2000 (Molecule &mol)
{
   MoleculeRGroups *rgroups = mol.getRGroups();

   if (rgroups == 0)
      rgroups = mol.createRGroups();

   // read groups
   while (!_scanner.isEOF())
   {
      char chars[5];

      chars[4] = 0;
      _scanner.readCharsFix(4, chars);

      if (strncmp(chars, "$RGP", 4) == 0)
      {
         _scanner.skipString();
         _scanner.skipSpace();

         int rgroup_idx = _scanner.readInt() - 1;
         RGroup &rgroup = rgroups->getRGroup(rgroup_idx);

         _scanner.skipString();
         while (!_scanner.isEOF())
         {
            char rgp_chars[6];
            rgp_chars[5] = 0;
            _scanner.readCharsFix(5, rgp_chars);

            if (strncmp(rgp_chars, "$CTAB", 5) == 0)
            {
               _scanner.skipString();
               AutoPtr<Molecule> fragment(new Molecule());

               MolfileLoader loader(_scanner);

               loader._readCtabHeader();
               fragment->createQueryAtoms();
               fragment->createQueryBonds();
               loader._readCtab2000(fragment.ref(), true);
               if (loader._rgfile)
                  loader._readRGroups2000(fragment.ref());
               loader._postLoad(fragment.ref(), true);

               rgroup.fragments.add(fragment.release());
            }
            else if (strncmp(rgp_chars, "$END ", 5) == 0)
            {
               rgp_chars[3] = 0;
               _scanner.readCharsFix(3, rgp_chars);

               _scanner.skipString();
               if (strncmp(rgp_chars, "RGP", 3) == 0)
                  break;
            }
            else
               _scanner.skipString();
         }
      }
      else if (strncmp(chars, "$END", 4) == 0)
      {
         chars[4] = 0;
         _scanner.readCharsFix(4, chars);
         _scanner.skipString();
         if (strncmp(chars, " MOL", 4) == 0)
            break;
      } else
         _scanner.skipString();
   }
}

void MolfileLoader::_readCtab3000 (Molecule &mol, bool query)
{
   QS_DEF(Array<char>, str);
   MoleculeRGroups *rgroups = mol.getRGroups();
   MoleculeRGroupFragment *rgroup_fragment = mol.getRGroupFragment();

   _scanner.readString(str, true);
   if (strncmp(str.ptr(), "M  V30 BEGIN CTAB", 17) != 0)
      throw Error("error reading CTAB block header");

   str.clear_resize(14);
   _scanner.read(14, str.ptr());
   if (strncmp(str.ptr(), "M  V30 COUNTS ", 14) != 0)
      throw Error("error reading COUNTS line");

   int i, nsgroups, n3d, chiral_int;

   _scanner.readString(str, true);
   if (sscanf(str.ptr(), "%d %d %d %d %d",
              &_atoms_num, &_bonds_num, &nsgroups, &n3d, &chiral_int) < 5)
      throw Error("error parsing COUNTS line");

   _chiral = (chiral_int != 0);

   _init();

   _scanner.readString(str, true);
   if (strncmp(str.ptr(), "M  V30 BEGIN ATOM", 14) != 0)
      throw Error("Error reading ATOM block header");

   for (i = 0; i < _atoms_num; i++)
   {
      _readMultiString(str);
      BufferScanner strscan(str.ptr());

      float x, y, z;
      int aamap;
      int iso = 0;
      int label = 0;
      int irflag = 0;
      int ecflag = 0;

      strscan.readInt1(); // atom index -- ignored

      QueryAtom query_atom_pre(0);
      QS_DEF(Array<char>, buf);

      strscan.readWord(buf, " [");

      char stopchar = strscan.readChar();

      if (stopchar == '[')
      {
         if (!query)
            throw Error("atom list is allowed only for queries");

         if (buf[0] == 0)
            query_atom_pre.type = QUERY_ATOM_LIST;
         else if (strcmp(buf.ptr(), "NOT") == 0)
            query_atom_pre.type = QUERY_ATOM_NOTLIST;
         else
            throw Error("bad word: %s", buf.ptr());

         bool was_a = false, was_q = false;

         while (1)
         {
            strscan.readWord(buf, ",]");
            stopchar = strscan.readChar();

            if (was_a)
               throw Error("'A' inside atom list, if present, must be single");
            if (was_q)
               throw Error("'Q' inside atom list, if present, must be single");

            if (buf.size() == 2 && buf[0] == 'A')
            {
               was_a = true;
               query_atom_pre.type = QUERY_ATOM_A;
            }
            else if (buf.size() == 2 && buf[0] == 'Q')
            {
               was_q = true;
               query_atom_pre.type = QUERY_ATOM_Q;
            }
            else
            {
               BufferScanner labelscan(buf.ptr());
               query_atom_pre.addToList(Element::read(labelscan));
            }

            if (stopchar == ']')
               break;
         }
      }
      else if (buf.size() == 2 && buf[0] == 'D')
      {
         label = ELEM_H;
         iso = 2;
      }
      else if (buf.size() == 2 && buf[0] == 'T')
      {
         label = ELEM_H;
         iso = 3;
      }
      else if (buf.size() == 2 && buf[0] == 'Q')
         query_atom_pre.type = QUERY_ATOM_Q;
      else if (buf.size() == 2 && buf[0] == 'A')
         query_atom_pre.type = QUERY_ATOM_A;
      else if (buf.size() == 2 && buf[0] == 'X' && !treat_x_as_pseudoatom)
         query_atom_pre.type = QUERY_ATOM_X;
      else if (buf.size() == 3 && buf[0] == 'R' && buf[1] == '#')
      {
         if (!query)
            throw Error("rgroup is allowed only for queries");

         if (mol.getRGroups() == 0)
            rgroups = mol.createRGroups();

         query_atom_pre.type = QUERY_ATOM_RGROUP;
         label = 0;
      }
      else
      {
         label = Element::fromString2(buf.ptr());
      }

      strscan.skipSpace();
      x = strscan.readFloat();
      strscan.skipSpace();
      y = strscan.readFloat();
      strscan.skipSpace();
      z = strscan.readFloat();
      strscan.skipSpace();
      aamap = strscan.readInt1();

      int idx = mol.addAtom(label);

      if (label == -1)
      {
         _preparePseudoAtomLabel(buf);
         mol.setPseudoAtom(idx, buf.ptr());
      }

      Atom &atom = mol.getAtom2(idx);
      QueryAtom *query_atom = 0;

      if (query)
      {
         query_atom = &mol.getQueryAtom2(idx);
         *query_atom = query_atom_pre;
      }
      else if (query_atom_pre.type != 0)
         throw Error("can't have query atoms in non-query molecule");

      atom.pos.set(x, y, z);
      atom.isotope = iso;

      if (reaction_atom_mapping != 0)
         reaction_atom_mapping->at(idx) = aamap;

      while (!strscan.isEOF())
      {
         QS_DEF(Array<char>, prop_arr);

         strscan.skipSpace();
         strscan.readWord(prop_arr, "=");

         strscan.skip(1);
         const char *prop = prop_arr.ptr();

         if (strcmp(prop, "CHG") == 0)
            atom.charge = strscan.readInt1();
         else if (strcmp(prop, "RAD") == 0)
            atom.radical = strscan.readInt1();
         else if (strcmp(prop, "CFG") == 0)
         {
            strscan.readInt1();
            //int cfg = strscan.readInt1();

            //if (cfg == 3)
            //   _stereocenter_types[idx] = MoleculeStereocenters::ATOM_ANY;
         }
         else if (strcmp(prop, "MASS") == 0)
            atom.isotope = strscan.readInt1();
         else if (strcmp(prop, "VAL") == 0)
         {
            atom.explicit_valence = true;
            atom.valence = strscan.readInt1();
         }
         else if (strcmp(prop, "HCOUNT") == 0)
         {
            if (query_atom != 0)
               query_atom->hydrogens = strscan.readInt1();
            else
               throw Error("nonzero atom hydrogens count is allowed only for queries");
         }
         else if (strcmp(prop, "STBOX") == 0)
            _stereo_care_atoms[idx] = strscan.readInt1();
         else if (strcmp(prop, "INVRET") == 0)
            irflag = strscan.readInt1();
         else if (strcmp(prop, "EXACHG") == 0)
            ecflag = strscan.readInt1();
         else if (strcmp(prop, "SUBST") == 0)
         {
            if (query_atom != 0)
               query_atom->substitution_count = strscan.readInt1();
            else
               throw Error("substitution count is allowed only for queries");
         }
         else if (strcmp(prop, "UNSAT") == 0)
         {
            if (query_atom != 0)
               query_atom->unsaturation = (strscan.readInt1() > 0);
            else
               throw Error("unsaturation flag is allowed only for queries");
         }
         else if (strcmp(prop, "RBCNT") == 0)
         {
            if (query_atom != 0)
               query_atom->ring_bond_count = strscan.readInt1();
            else
               throw Error("ring bond count is allowed only for queries");
         }
         else if (strcmp(prop, "RGROUPS") == 0)
         {
            if (query_atom == 0)
               throw Error("rgroup is allowed only for queries");

            if (mol.getRGroups() == 0)
               rgroups = mol.createRGroups();

            int n_rg, rg_idx;

            strscan.skip(1); // skip '('
            n_rg = strscan.readInt1();
            while (n_rg-- > 0)
            {
               rg_idx = strscan.readInt1() - 1;
               rgroups->addRGroupOnSite(idx, rg_idx);
            }
         }
         else if (strcmp(prop, "ATTCHPT") == 0)
         {
            if (!query)
               throw Error("rgroup attachment points are allowed only for queries");

            if (mol.getRGroupFragment() == 0)
               rgroup_fragment = mol.createRGroupFragment();

            int att_type = strscan.readInt1();

            if (att_type == -1)
               att_type = 3;

            for (int att_idx = 0; (1 << att_idx) <= att_type; att_idx++)
               if (att_type & (1 << att_idx))
                  rgroup_fragment->addAttachmentPoint(att_idx, idx);
         }
         else if (strcmp(prop, "ATTCHORD") == 0)
         {
            if (query_atom != 0)
            {
               if (mol.getRGroups() == 0)
                  rgroups = mol.createRGroups();

               int n_items, nei_idx, att_type;

               strscan.skip(1); // skip '('
               n_items = strscan.readInt1();
               while (n_items-- > 0)
               {
                  nei_idx = strscan.readInt1();
                  att_type = strscan.readInt1();
                  rgroups->setAttachmentOrder(idx, att_type - 1, nei_idx);
               }
            }
            else
               throw Error("rgroup is allowed only for queries");
         }
         else
         {
            throw Error("unsupported property of CTAB3000: %s", prop);
         }
      }

      if (reaction_atom_inversion != 0)
         reaction_atom_inversion->at(idx) = irflag;

      if (reaction_atom_exact_change != 0)
         reaction_atom_exact_change->at(idx) = ecflag;
   }

   _scanner.readString(str, true);
   if (strncmp(str.ptr(), "M  V30 END ATOM", 15) != 0)
      throw Error("Error reading ATOM block footer");

   _scanner.readString(str, true);
   if (strncmp(str.ptr(), "M  V30 BEGIN BOND", 17) != 0)
   {  
      if (_bonds_num > 0)
         throw Error("Error reading BOND block header");
   }
   else
   {
      for (i = 0; i < _bonds_num; i++)
      {
         int reacting_center = 0;

         _readMultiString(str);
         BufferScanner strscan(str.ptr());

         strscan.readInt1(); // bond index -- ignored

         int type = strscan.readInt1();
         int beg = strscan.readInt1() - 1;
         int end = strscan.readInt1() - 1;
         int query_type = 0;

         switch (type)
         {
            case BOND_SINGLE: 
            case BOND_DOUBLE:
            case BOND_TRIPLE:
            case BOND_AROMATIC:
               break;
            case QUERY_BOND_SINGLE_OR_DOUBLE:
               if (!query)
                  throw Error("'single or double' bonds are allowed only for queries");
               query_type = type;
               type = 0;
               break;
            case QUERY_BOND_SINGLE_OR_AROMATIC:
               if (!query)
                  throw Error("'single or aromatic' bonds are allowed only for queries");
               query_type = type;
               type = 0;
               break;
            case QUERY_BOND_DOUBLE_OR_AROMATIC:
               if (!query)
                  throw Error("'double or aromatic' bonds are allowed only for queries");
               query_type = type;
               type = 0;
               break;
            case QUERY_BOND_ANY:
               if (!query)
                  throw Error("'any' bonds are allowed only for queries");
               query_type = type;
               type = 0;
               break;
            default:
               throw Error("unknown bond type: %d", type);
         }

         int idx = mol.addBond(beg, end, type);
         QueryBond *query_bond = 0;

         if (query)
         {
            query_bond = &mol.getQueryBond2(idx);
            query_bond->type = query_type;
         }

         while (!strscan.isEOF())
         {
            QS_DEF(Array<char>, prop);

            strscan.readWord(prop, "=");
            strscan.skip(1);

            if (strcmp(prop.ptr(), "CFG") == 0)
            {
               int n = strscan.readInt1();

               if (n == 1)
                  _bond_directions[i] = MoleculeStereocenters::BOND_UP;
               else if (n == 3)
                  _bond_directions[i] = MoleculeStereocenters::BOND_DOWN;
               else if (n == 2)
               {
                  _bond_directions[i] = MoleculeStereocenters::BOND_EITHER;
                  _stereocenter_types[beg] = MoleculeStereocenters::ATOM_ANY;
               }
               else
                  throw Error("unknown bond CFG=%d", n);

               if (_stereocenter_types[beg] == 0)
               {
                  _stereocenter_types[beg] = _chiral ? MoleculeStereocenters::ATOM_ABS : MoleculeStereocenters::ATOM_AND;
                  _stereocenter_groups[beg] = 1;
               }
            }
            else if (strcmp(prop.ptr(), "STBOX") == 0)
            {
               if (query_bond != 0)
                  query_bond->stereo_care = 1;
               else
                  throw Error("stereo care box is allowed only for queries");
               query_bond->stereo_care = (strscan.readInt1() != 0);
            }
            else if (strcmp(prop.ptr(), "TOPO") == 0)
            {
               if (query_bond != 0)
                  query_bond->topology = strscan.readInt1();
               else
                  throw Error("bond topology setting is allowed only for queries");
            }
            else if (strcmp(prop.ptr(), "RXCTR") == 0)
               reacting_center = strscan.readInt1();
         }

         if (reaction_bond_reacting_center != 0)
            reaction_bond_reacting_center->at(idx) = reacting_center;
      }

      _scanner.readString(str, true);
      if (strncmp(str.ptr(), "M  V30 END BOND", 15) != 0)
         throw Error("Error reading BOND block footer");

      _scanner.readString(str, true);
   }

   if (highlighting != 0)
      highlighting->init(mol);

   if (strncmp(str.ptr(), "M  V30 BEGIN COLLECTION", 23) == 0)
   {
      while (1)
      {
         _readMultiString(str);

         if (strncmp(str.ptr(), "END COLLECTION", 14) == 0)
            break;

         BufferScanner strscan(str.ptr());
         char coll[14];

         strscan.readCharsFix(13, coll);
         coll[13] = 0;

         int stereo_type = 0;
         int stereo_group = 0;
         int n = 0;

         if (strcmp(coll, "MDLV30/STERAC") == 0)
            stereo_type = MoleculeStereocenters::ATOM_AND;
         else if (strcmp(coll, "MDLV30/STEREL") == 0)
            stereo_type = MoleculeStereocenters::ATOM_OR;
         else if (strcmp(coll, "MDLV30/STEABS") == 0)
            stereo_type = MoleculeStereocenters::ATOM_ABS;
         else if (strcmp(coll, "MDLV30/HILITE") == 0)
         {
            if (highlighting != 0)
            {
               QS_DEF(Array<char>, what);

               strscan.skipSpace();
               strscan.readWord(what, " =");

               if (strcmp(what.ptr(), "ATOMS") == 0)
               {
                  strscan.skip(2); // =(
                  n = strscan.readInt1();
                  while (n -- > 0)
                     highlighting->onVertex(strscan.readInt1() - 1);
               }
               else if (strcmp(what.ptr(), "BONDS") == 0)
               {
                  strscan.skip(2); // =(
                  n = strscan.readInt1();
                  while (n -- > 0)
                     highlighting->onEdge(strscan.readInt1() - 1);
               }
               else
                  throw Error("unknown hightlighted object: %s", what.ptr());
            }
            
            continue;
         }
         else
            throw Error("unknown collection: %s", coll);

         if (stereo_type == MoleculeStereocenters::ATOM_OR ||
             stereo_type == MoleculeStereocenters::ATOM_AND)
            stereo_group = strscan.readInt1();
         else
            strscan.skip(1);

         strscan.skip(7); // ATOMS=(
         n = strscan.readInt1();
         while (n -- > 0)
         {
            int atom_idx = strscan.readInt1() - 1;
            
            _stereocenter_types[atom_idx] = stereo_type;
            _stereocenter_groups[atom_idx] = stereo_group;
         }
      }
      _scanner.readString(str, true);
   }

   if (strncmp(str.ptr(), "M  V30 END CTAB", 15) != 0)
      throw Error("error reading CTAB block footer");

}

void MolfileLoader::_preparePseudoAtomLabel (Array<char> &pseudo)
{
   // if the string is quoted, unquote it
   if (pseudo.size() > 2 && pseudo[0] == '\'' && pseudo[pseudo.size() - 2] == '\'')
   {
      pseudo.remove(pseudo.size() - 2);
      pseudo.remove(0);
   }

   // make the first letter capital and the others lowercase
   //pseudo.tolower();
   //if (pseudo.size() > 0 && isalpha(pseudo[0]))
   //   pseudo[0] = toupper(pseudo[0]);

   if (pseudo.size() <= 1)
      throw Error("empty pseudo-atom");
}

void MolfileLoader::_readMultiString (Array<char> &str)
{
   QS_DEF(Array<char>, tmp);

   str.clear();
   tmp.clear_resize(7);

   while (1)
   {
      bool to_next = false;

      _scanner.read(7, tmp.ptr());
      if (strncmp(tmp.ptr(), "M  V30 ", 7) != 0)
         throw Error("error reading multi-string in CTAB v3000");

      _scanner.readString(tmp, true);

      if (tmp[tmp.size() - 2] == '-')
      {
         tmp[tmp.size() - 2] = 0;
         tmp.pop();
         to_next = true;
      }
      str.appendString(tmp.ptr(), true);
      if (!to_next)
         break;
   }
}

void MolfileLoader::_readRGroups3000 (Molecule &query)
{
   QS_DEF(Array<char>, str);

   MoleculeRGroups *rgroups = query.getRGroups();

   if (rgroups == 0)
      rgroups = query.createRGroups();

   while (!_scanner.isEOF())
   {
      _scanner.readString(str, true);

      if (strncmp(str.ptr(), "M  V30 BEGIN RGROUP", 19) == 0)
      {
         _rgfile = true;

         int rg_idx;

         if (sscanf(str.ptr(), "M  V30 BEGIN RGROUP %d", &rg_idx) != 1)
            throw Error("can not read rgroup index");

         RGroup &rgroup = rgroups->getRGroup(rg_idx - 1);

         _readMultiString(str);

         BufferScanner strscan(str.ptr());

         if (strncmp(str.ptr(), "RLOGIC", 6) != 0)
            throw Error("Error reading RGROUP block");

         strscan.skip(7);
         rgroup.if_then = strscan.readInt1() - 1;
         rgroup.rest_h = strscan.readInt1();

         QS_DEF(Array<char>, occ);

         strscan.readString(occ, true);
         _readRGroupOccurenceRanges(occ.ptr(), rgroup.occurence);

         while (!_scanner.isEOF())
         {
            int pos = _scanner.tell();

            _scanner.readString(str, true);
            if (strcmp(str.ptr(), "M  V30 BEGIN CTAB") == 0)
            {
               _scanner.seek(pos, SEEK_SET);
               AutoPtr<Molecule> fragment(new Molecule());
               fragment->createQueryAtoms();
               fragment->createQueryBonds();
               MolfileLoader loader(_scanner);
               loader._readCtab3000(fragment.ref(), true);
               loader._postLoad(fragment.ref(), true);
               rgroup.fragments.add(fragment.release());
            }
            else if (strcmp(str.ptr(), "M  V30 END RGROUP") == 0)
               break;
            else
               throw Error("unexpected string in rgroup: %s", str.ptr());
         }

      }
      else if (strncmp(str.ptr(), "M  END", 6) == 0)
         break;
      else
         throw Error("unexpected string in rgroup: %s", str.ptr());
   }
}

void MolfileLoader::_init ()
{
   _stereo_care_atoms.clear_resize(_atoms_num);
   _stereo_care_atoms.zerofill();
   _stereocenter_types.clear_resize(_atoms_num);
   _stereocenter_types.zerofill();
   _stereocenter_groups.clear_resize(_atoms_num);
   _stereocenter_groups.zerofill();
   _bond_directions.clear_resize(_bonds_num);
   _bond_directions.zerofill();
   _ignore_cistrans.clear_resize(_bonds_num);
   _ignore_cistrans.zerofill();

   if (reaction_atom_mapping != 0)
      reaction_atom_mapping->clear_resize(_atoms_num);
   if (reaction_atom_inversion != 0)
      reaction_atom_inversion->clear_resize(_atoms_num);
   if (reaction_atom_exact_change != 0)
      reaction_atom_exact_change->clear_resize(_atoms_num);
   if (reaction_bond_reacting_center != 0)
      reaction_bond_reacting_center->clear_resize(_bonds_num);

   if (highlighting != 0)
      highlighting->clear();
}

void MolfileLoader::_createStereocenters (Molecule &mol)
{
   int i;

   for (i = mol.vertexBegin(); i != mol.vertexEnd(); i = mol.vertexNext(i))
      if (_stereocenter_types[i] > 0)
         break;

   if (i == mol.vertexEnd())
      return;

   MoleculeStereocenters & stereocenters = mol.getStereocenters();

   stereocenters.buildFromBonds(_stereocenter_types.ptr(),
      _stereocenter_groups.ptr(), _bond_directions.ptr(), ignore_stereocenter_errors);
}
